Body mass dynamics of breeding male pectoral sandpipers

Aims and sample overview

We want to understand how body mass varies within and among male pectoral sandpipers and how these dynamics associate with reproductive success and accelerometry. Our dataset utilizes five seasons of repeated measures of breeding males collected between 2005 to 2009 at the Barrow study site in arctic Alaska, totaling 151 repeated measures from 71 males. During the 2005 to 2009 period, we conducted intense monitoring of the breeding activity, mating dynamics, and paternity outcomes of males - enableing us to assess how body mass changes relate to reproductive effort and success. Furthermore, in 2023, we conducted a concerted effort to recapture males that were deployed with GPS/Accelerometer tags, totaling 91 repeated measures from 41 males. This latter dataset lacks reproductive information, but the tags allow us to assess how individual energy expenditure (measured through accelerometry) is linked to within-individual body mass variation.

Data wrangling

Standardize captures by year

Despite the short season, early arriving males likely face different conditions than later arriving males that may impact body mass dynamics during their tenure at Barrow. For example, early arriving males may face different levels of intra-sexual competition, food availability, or ambient climate-conditions than later arriving males which would influence metabolic demands, and hence shifts in repeated measures of body mass. To standardize measurements to the year-specific phenology of males, we first extract the first captures of all males in years 2005 to 2009 from the PESAatBARROW.CAPTURES database (i.e., even those that were captured only once). Then, for each capture, we calculate the time since 00:00:00 on January 1 of a capture’s given year to acquire a Julian-esque timestamp.

Then we standardize capture dates by the mean capture date for each year

Evaluate static morphometric measures of body structure with a PCA

Body mass changes may differ according to the structural size of an individual and so in our model we will need to control for structural size. In the PESAatBARROW.CAPTURES database we have 4 structural size measures collected for each individual: culmen, totalHead, tarsus, and wing. Here we run a Principle Component Analysis of these 4 measures to assess their collinearity, then we evaluate how each measure correlates with body mass compared to PC1 of the PCA. The results of the PCA show that PC1 captures ~81% of the variance in the 4 measures, however when we compare the correlations of PC1 with weight vs. wing with weight, we can see that we do not gain much more of inference with PC1. Thus we proceed with using wing as our structural size control as it is more parsimonious.

Importance of components:
                          Comp.1    Comp.2     Comp.3     Comp.4
Standard deviation     3.7724066 1.5484817 0.88062916 0.46171432
Proportion of Variance 0.8077777 0.1361028 0.04401908 0.01210045
Cumulative Proportion  0.8077777 0.9438805 0.98789955 1.00000000

Extract tenures

visualize annual variation in tenure data for the full and reduced datasets

Within-individual centering: Disentagling within- vs. between individual changes in body mass over time

When assessing how population values of body mass change over time or differ between individuals, it is important to acknowledge that both within- and between-individual processes might be underlying causal mechanisms (van de Pol & Verhulst, 2006). We therefore need to decompose the within- vs. between-individual components of body mass change in the male PESA population. Population-level correlations might result from seasonal variation in the arrival (or departure) of males of differing body mass. This is known as ‘selective appearance’ (or ‘selective disappearance’, in the case of last measures). In our study of PESA body mass, one could hypothesize that there is selective appearance of large males: heavier males arrive earlier than lighter males, due to faster spring migration, etc. Moreover, one could also hypothesize that there is selective disappearance of light males: lighter males may depart Barrrow sooner than heavy males due to their disadvantaged local intra-sexual competitive abilities. These effects result in correlations between tenure and body mass in our cross-sectional analysis of body mass, without reflecting within-individual changes. Within-individual changes in body mass can be caused by factors intrinsic to the individual, such as physical state, reproductive tactics, or by extrinsic factors, such as changes in the availability of food or fertile females on the study site. Within- and between-individual changes are not mutually exclusive: they may go in the same direction (Fig. 1A), or they may go in opposite directions with the between-individual effect masking within-individual patterns of an increase (Fig. 1B) or a decline (Fig. 1C) in body mass.

Figure 1. Schematic showing a hypothetical population of three individuals (faint lines) of different body mass dynamics, compared to ordinary regression lines fitted through all individuals (bold-dashed lines). A few possible scenarios shown here are (A) within-individual increase in body mass and selective appearance of heavier individuals over the season, (B) within-individual increase in body mass and selective disappearance of heavy individuals over the season, and (C) within-individual decline in body mass and selective disappearance of lighter males over the season. Modified from van de Pol & Verhulst, 2006, American Naturalist 167:766-773.

To control for selective appearance and disappearance of males differing in body mass we fit our mixed-effects model with ‘date first measured’ and ‘date last measured’ as fixed effects – a method that estimates between-individual seasonal effects introduced by selective disappearance and appearance. The raw distributions of ‘date first measured’ and ‘date last measured’ are shown here:

We model within-individual temporal changes in body mass by fitting a within-individual deviation score for date (henceforth ‘date-deviance’), calculated for individual \(i\) at date \(j\) as: \(date_{ij} - {Date FirstMeasured}_i\) (van de Pol & Verhulst, 2006a; Snijders & Bosker, 2011) – essentially describing the amount of time since the bird’s first encounter in the population (assumed to be a good proxy for the amount of time in Barrow since arrival).

Mixed effects model of within- vs. between-individual body mass dynamics

Our model evaluates repeated measures of body mass (\(r_{ij}\)) taken from individual \(i\) at time \(j\). The random intercept term \(u_{0i}\) and the residual error term \(e_{0ij}\) are assumed to be drawn from a Gaussian distribution with mean 0 and variance \(\sigma^2_u\) and \(\sigma^2_e\), respectively. Our model specifically tests for within-individual change (\(\beta_W \times DateDeviance_{ij}\)) in the presence of a selective appearance effect (\(\beta_F \times FirstDate_i\)) and a selective disappearance effect (\(\beta_L \times LastDate_i\)). Thus the effect of date on first measure (\(\beta_F\)) estimates the additional effect to selective appearance given the estimated within-individual change in body mass with time. Wing (\(\beta_G \times Wing_{i}\)) is included in our model to control for the structural size of individuals and year and individual identity were included as a random effects: \[ r_{ij} = \beta_0 + \beta_G \times Wing_{i} + \beta_W \times DateDeviance_{ij} + \beta_F \times FirstDate_i + \beta_L \times LastDate_i + u_{0i} + e_{0ij} \]

Code

mod_weight <-
  lmer(weight ~ wing + date_deviance + first_date + last_date +
        (1 | capture_id) + (1 | year_),
       data = cap_05_09_23_std_pca_ten)

To evaluate uncertainty in our parameter estimates we simulated 1000 parametric bootstraps via “broom.mixed::tidy” and the “partR2::partR2” function (Stoffel et al., 2020). Likewise we derived individual-level repeatabilities (i.e., intra-class correlation) by simulating 1000 parametric bootstraps of the mixed-effect model using “rptR::rpt”. We report fixed effects as standardized regression coefficients (i.e., beta weights) and repeatability as the ‘adjusted repeatability’ – interpreted as the repeatability of a given hierarchical group after controlling for fixed effects (Nakagawa & Schielzeth, 2010).

First let’s plot the residuals against the fitted values. Looks good.

Then let’s explore the effect sizes and predictions.

wing and weight are tightly correlated as expected
strong negative within-individual effect of season on body mass (after controlling for wing length)
strong negative between-individual effect of season on body mass (after controlling for wing length): “selective appearance”
no between-individual effect of season on body mass (after controlling for wing length): no “selective disappearance”

male Pectoral Sandpiper body mass dynamics	Mean estimate	95% confidence interval
Fixed effects 𝛽 (standardized)
Wing length	0.41	[0.29, 0.52]
Within ind. temporal change	−0.21	[-0.31, -0.11]
Between ind. effect of season (first measure)	−0.23	[-0.38, -0.06]
Between ind. effect of season (last measure)	0.03	[-0.13, 0.19]
Partitioned 𝑹²
Total Marginal 𝑹²	0.23	[0.15, 0.34]
Wing length	0.16	[0.08, 0.28]
Within ind. temporal change	0.03	[0, 0.16]
Between ind. effect of season (first measure)	0.00	[0, 0.13]
Between ind. effect of season (last measure)	0.00	[0, 0.13]
Total Conditional 𝑹²	0.51	[0.39, 0.64]
Random effects 𝜎²
Individual	13.77	[5.42, 23.32]
Year	4.20	[0, 15.04]
Residual	30.68	[24.15, 39.48]
Adjusted repeatability 𝑟
Individual	0.28	[0.13, 0.44]
Year	0.08	[0, 0.23]
Residual	0.64	[0.48, 0.8]
Sample sizes 𝑛
Years	6
Individuals	113
Observations	242

Relationship between body mass and wing length:

as expected, strong positive relationship between body mass and wing length

Within-individual change in body mass over season

decline over season: males tend to loose 1g every ~3.8 days

Inspecting the increases and decreases

Inspecting the males with > 2 measures

Between-individual effect of first date (i.e., selective appearance)

strong negative relationship: heavier males arrive earlier than lighter males

Between-individual effect of last date (i.e., selective disappearance)

no effect of last measure date: departure from Barrow is not related to body mass

Linking body mass dynamics to paternity

The next step is to link the body mass dynamics to fitness traits of interest, specifically the number of female mates acquired, the the number of young sired, and a binary form of paternity (i.e., sired at least one offspring vs. none).

Annual variation in polygyny data for repeated measures dataset

Annual variation in paternity data for repeated measures dataset

Annual variation in paternity data (binary) for repeated measures dataset

Next we need to transform our repeated measures of body mass into measure of change (e.g., \(\Delta BodyMass\)), which is essentially calculating the difference between the first and last measure. We will also need to bind the paternity data to the dataframe. We are left with a single row per male that allows for models assessing how the change in body condition is associated with polygyny and paternity.

check how delta weight varies according to the time difference between measures: no relationship

Code

mod_abs_delta_weight_date_deviance <- 
  lmer(abs_delta_weight ~ date_deviance + (1 | year_),
        data = cap_05_09_std_pca_ten_pat_delta_res)

Mixed effects models of mating dynamics, tenure, and body mass

Is change in body mass associated with polygyny?

Heavy males have more females, as do those that have longer tenure. No effect of weight change on polygyny.

Code

mod_n_females_weight1_delta_tenure <-
  glmer(N_females_0 ~ weight_1 + delta_weight + tenure + (1 | year_),
        data = cap_05_09_std_pca_ten_pat_delta_res, family = poisson)

Characteristic	log(IRR)¹	95% CI¹	p-value
(Intercept)	-7.7	-12, -3.3	<0.001
weight_1	0.07	0.02, 0.11	0.004
delta_weight	0.02	-0.02, 0.07	0.3
tenure	0.05	0.01, 0.08	0.011
¹ IRR = Incidence Rate Ratio, CI = Confidence Interval

Is change in body mass associated with paternity?

Heavy males sire more young, as do those that have longer tenure. Males that minimize weight loss have higher paternity.

Code

mod_n_young_weight1_delta_tenure <-
  glmer(N_young_0 ~ weight_1 + delta_weight + tenure + (1 | year_),
        data = cap_05_09_std_pca_ten_pat_delta_res, family = poisson)

Characteristic	log(IRR)¹	95% CI¹	p-value
(Intercept)	-7.5	-10, -5.1	<0.001
weight_1	0.07	0.05, 0.10	<0.001
delta_weight	0.03	0.01, 0.06	0.006
tenure	0.05	0.03, 0.07	<0.001
¹ IRR = Incidence Rate Ratio, CI = Confidence Interval

Is change in body mass associated with paternity probability?

Males with longer tenure have the highest chance to sire at least one young. No effect of initial wieght or body mass change on paternity probability.

Code

mod_weight1_delta_tenure_pat_bin <-
  glmer(cbind(pat, no_pat) ~ weight_1 + delta_weight + tenure + (1 | year_),
        data = cap_05_09_std_pca_ten_pat_delta_res, family = "binomial")

Characteristic	log(OR)¹	95% CI¹	p-value
(Intercept)	-8.2	-17, 0.09	0.053
weight_1	0.06	-0.02, 0.15	0.13
delta_weight	0.00	-0.07, 0.07	>0.9
tenure	0.09	0.01, 0.16	0.023
¹ OR = Odds Ratio, CI = Confidence Interval

Linking Overall Body Dynamic Acceleration to change in body mass

 [1] "049C" "0453" "0447" "048D" "048E" "0468" "044D" "0469" "0437" "043F"
[11] "0433" "0438" "047E" "0443" "0484" "043f" "0457" "0477" "0467" "048F"
[21] "044A" "0449" "044F" "0473" "0472" "047B" "0475" "0462" "0494" "0432"
[31] "045B" "0463" "045C" "0493" "0445" "049B" "0444" "0485" "0455" "0434"
[41] "045D" "045A" "0461" "048C" "0492" "047F" "0460" "046E" "0451" "0490"
[51] "0459" "0495" "0458" "046D" "0454" "03E3" "0450" "043E" "043B" "049E"
[61] "047A" "0465" "0442" "0440" "0478" "0435" "046B" "0486" "0488" "0479"
[71] "046A" "046C" "043A" "0491" "044E" "044B" "049D" "0487" "0476" "0471"
[81] "0466" "044C" "0481" "0441" "0452" "0482" "0448" "0456"

# A tibble: 1 × 5
  capture_id tagID start_tag_data      end_tag_data         dead
  <chr>      <chr> <dttm>              <dttm>              <dbl>
1 69407_2023 048A  2023-06-12 14:20:00 2023-06-17 01:52:00     1

         date cap_start   caught released form_id author gps_id gps_point    ID
1: 2023-06-16  15:40:00 15:52:00       NA     103     BK      1        59 69407
   UL UR   LL       LR tagID tag_action recapture tarsus culmen total_head wing
1:  M  W <NA> DB,PI,PI  <NA>        off         1     NA     NA         NA   NA
   weight wt_w_tag swab dead camID photoID    comments nov  pk
1:  98.26        0    1    1    NA      NA brain bird.   0 132

inspect temporal changes in ODBA for the males that we have repeated measures of body mass (windows of ODBA reflect the time between first and last capture). Black points mark body mass measuring events.

summarise ODBA data for each individual between body mass measures and convert delta weight to a daily rate to account for varied time windows between body mass measures).

Average ODBA model

model assessing the relationship between average ODBA and change in body mass (converted to daily rate to account for varied time windows between body mass measures).

Code

# non-standardized model
mod_avgODBA_delta <-
  lm(daily_delta_weight ~ avg_ODBA + weight_1 + first_date,
        data = cap_05_09_23_std_pca_delta_odba)

# standardized model
mod_avgODBA_delta_ <-
  lm(daily_delta_weight ~ avg_ODBA_ + weight_1_ + first_date_,
        data = cap_05_09_23_std_pca_delta_odba)

check for colinearity among predictor variables (looks good!)

strong positive effect of ODBA on change in mass (more active males experience less mass change over time)
strong negative effect of initial mass on change in mass (heavier males experience more mass change over time)
weak negative effect of first capture date on change in mass (late arriving males experience more mass change over time) Note: table shows standardized effect sizes.

Characteristic	Beta	95% CI¹	p-value
(Intercept)	-0.69	-1.0, -0.34	<0.001
avg_ODBA_	0.46	0.10, 0.83	0.014
weight_1_	-0.99	-1.4, -0.63	<0.001
first_date_	-0.33	-0.70, 0.03	0.071
¹ CI = Confidence Interval

Bootstrap of the average ODBA model

Average ODBA model (percent change in body mass)

model assessing the relationship between average ODBA and percent change in body mass (converted to daily rate to account for varied time windows between body mass measures).

Code

# non-standardized model
mod_avgODBA_delta_p <-
  lm(daily_delta_weight_perc ~ avg_ODBA + first_date + weight_1,
        data = cap_05_09_23_std_pca_delta_odba)

# standardized model
mod_avgODBA_delta_p_ <-
  lm(daily_delta_weight_perc ~ avg_ODBA_ + first_date_ + weight_1,
        data = cap_05_09_23_std_pca_delta_odba)

strong positive effect of ODBA on change in mass (more active males experience less mass change over time)
strong negative effect of initial mass on change in mass (heavier males experience more mass change over time)
weak negative effect of first capture date on change in mass (late arriving males experience more mass change over time) Note: table shows standardized effect sizes.

Characteristic	Beta	95% CI¹	p-value
(Intercept)	13	7.5, 18	<0.001
avg_ODBA_	0.47	0.12, 0.82	0.010
first_date_	-0.32	-0.67, 0.03	0.073
weight_1	-0.13	-0.18, -0.08	<0.001
¹ CI = Confidence Interval

Bootstrap of the average ODBA model

Geometric mean ODBA model

model assessing the relationship between geometric mean ODBA and change in body mass (converted to daily rate to account for varied time windows between body mass measures).

Code

# non-standardized model
mod_gmODBA_delta <-
  lm(daily_delta_weight ~ gm_ODBA + weight_1 + first_date,
        data = cap_05_09_23_std_pca_delta_odba)

# standardized model
mod_gmODBA_delta_ <-
  lm(daily_delta_weight ~ gm_ODBA_ + weight_1_ + first_date_,
        data = cap_05_09_23_std_pca_delta_odba)

strong positive effect of ODBA on change in mass (more active males experience less mass change over time)
strong negative effect of initial mass on change in mass (heavier males experience more mass change over time)
weak negative effect of first capture date on change in mass (late arriving males experience more mass change over time) Note: table shows standardized effect sizes.

Characteristic	Beta	95% CI¹	p-value
(Intercept)	-0.69	-1.1, -0.34	<0.001
gm_ODBA_	0.44	0.07, 0.81	0.020
weight_1_	-0.95	-1.3, -0.58	<0.001
first_date_	-0.35	-0.72, 0.01	0.056
¹ CI = Confidence Interval

Bootstrap of the geometric mean ODBA model